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QN^ We develop the hybrid Monte Carlo method for simulations of single off-lattice polymer 

chains. We discuss implementation and choice of simulation parameters in some detail. 
The performance of the algorithm is tested on models for homopolymers with short- or 
long-range self-repulsion, using chains with 16 < iV < 512 monomers. Without excessive 
<D fine tuning, we find that the computational cost grows as N^+^ with 0.64 < z' < 0.84. In 

addition, we report results for the scaling of the end-to-end distance, r"i^ ~ iV'^(lniV)~". 
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1 INTRODUCTION 



Monte Carlo methods is a well-established tool in the study of polymer models and many 
efficient algorithms have been developed [1]. To facilitate the study of large systems it is, 
nevertheless, of interest to continue the search for improved methods. In this paper we 
investigate the hybrid Monte Carlo method (HMC) [2], which can be used to simulate off- 
lattice chains. Although based on molecular dynamics, the method has, to our knowledge, 
not been applied to polymer models before. 

In HMC the system is evolved by using molecular dynamics in a fictitous time, which 
is supplemented with refreshments of the "momenta". The algorithm is made exact by 
incorporating a Metropolis type accept-reject step, which eliminates finite-step-size errors 
arising from the discretization of the equations of motion. The method has become widely 
used for simulations of lattice QCD with dynamical fermions. Important there is that a 
HMC update of the whole system requires only 0(1) calculations of the highly non-local 
Boltzmann weight. In this paper we apply HMC to single chains with interactions between 
all monomer pairs. A whole chain of this kind can be updated in a computer time of order 
N"^ . As an example, let us compare this with the behaviour of the popular pivot algorithm 
[3]. This algorithm was thoroughly tested for self-avoiding walks on a lattice in ref. [4], 
and was found very powerful for generating independent measurements of global quantities. 
However, the computer time required to generate one independent measuremnet of local 
quantities grows as or faster for the chains considered here, since each elementary pivot 
move scales as N"^ . The hope is that HMC can be used to improve on this, without losing 
good efficiency for global quantities. 

There are different ways to implement HMC, and the resulting efficiency varies. In our 
calculations we employ the so-called Fourier acceleration technique [5]. For linear chains 
with free endpoints, this is equivalent to performing ordinary HMC updates of the bond 
variables. We have tested this algorithm on four models with self-repulsion. We find that 
the computational cost of generating one independent measurement grows as iV2+^' with z' 
beween 0.64 and 0.84, and that it is similar for global and local quantities. In particular, 
these results suggest that thermalization can be carried out faster with HMC than with the 
pivot algorithm for large N . This property is important since thermalization can be the 
dominating cost with the pivot algorithm. 

The plan of this paper is as follows. In section 2 we present the models considered. In section 
3 we describe the algorithm and the tuning of simulation parameters. The measurements 
of the efficiency of the algorithm are discussed in section 4, where we also present results 
for the scaling of the end-to-end distance. A summary is given in section 5. 

2 THE MODELS 

Our tests of the performance of HMC have been carried out on models with long- or short- 
range self-repulsion. We first consider a polyelectrolyte in a solution. Here the repulsive 
interaction is taken to be of screened Coulomb type, and the presence of a salt concentration 
is parametrized through a finite Debye screening length. The full potential energy considered 
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where q is the monomer charge, e^eo the dielectric permittivity of the medium and 
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is the distance between monomers i and The inverse screening length k at salt concen 
tration is given by 

where is Avogadro's number and T the temperature. For convenience, all lengths will 
be measured in units of r"o = ( ^^^ eofe )^^^" "^^^ partition function then becomes 
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where k = tqk and T = [krQ)~^kBT . The endpoints of the chains are free and the S function 
is needed because of translational invariance. Our calculations have been performed using 
Cg = (unscreened Coulomb potential) or IM. The other parameters were taken to be 
ro = 6A, q = e, e = 78.3 and and T = 298K. This corresponds to T = 0.838 and to 
K = 1.992 for Cs = IM. These two sets of parameters were chosen to make a comparison 
possible with the recent results of refs. [6, 7]. These authors investigated the model using 
a Gaussian variational method, and performed also Monte Carlo calculations for N up to 
2048. 

The scaling with N of the end-to-end distance defines the critical index v through 

{riN) , N ^ oo . (5) 

Logarithmic corrections to this asymptotic relation may appear. The unscreened and 
screened Coulomb models represent different values oil/. The unscreened potential gives rise 
to rod-like behaviour with v = 1, while the Flory result v = 3/{2 + d) = 0.6 is approximately 
valid for the short-range screened potential. 

In addition to these two cases, we also consider two potentials of intermediate range. Here 
the repulsion energy decays as with A = 2 and 2.5, respectively, and the partition 
function is 
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where we put the dimensionless temperature parameter T = 1. This model was investigated 
in ref. [8] for 2 < A < 3. Using a Gaussian variational approach, these authors predicted 
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that V = 2/A, which differs from the Flory-type result v = 3/(2 + A). For A = 2, they 
further predicted the appearance of logarithmic corrections to eq. 5 of the form 

(r-iiv) ~iV^(lniV)-« (7) 

with = 2/A = 1 and a = 1/2. These results were tested numerically in ref. [9], using 
chains of size N < 120. The Monte Carlo data gave support to the result v = 2/X and 
to the presence of logarithmic corrections for A = 2. However, the precise value of the 
exponent a could not be determined from these data. In section 4, we will therefore use 
our results to get an improved estimate of a. We will find that the prediction of ref. [8] is 
in nice agreement with the data. 

Finally, we mention the virial theorem which provides a useful check of the Monte Carlo 
procedure, as discussed in refs. [6, 7]. For the partition functions in eqs. 4 and 6 this exact 
identity takes the forms 

iV-l -Rr- 

E V--^ ^ e--^-.) = 3T(iV-l) (8) 

i=l l<i<j<N *J l<i<j<N 

and 

(E^'i+i-^- E 4r) = 3T(iV - 1) , (9) 

i=l l<i<j<N ^ij 

respectively. 

3 THE ALGORITHM 

The HMC algorithm [2] is a general method for simulation of statistical systems with con- 
tinuous degrees of freedom and can be directly applied to the models considered here. Naive 
HMC updates of the monomer coordinates leads, however, to poor performance due to the 
slow evolution of long-wavelength modes. For linear chains with free endpoints there are 
two equivalent ways to alleviate this problem. One is, as in the usual Metropolis algorithm, 
to instead update the bond variables bi^ = Xi^i^ — Xi^. The other is to employ the Fourier 
acceleration technique [5]. 

3.1 HMC 

We begin with a brief description of ordinary HMC using the bond variables bi^. To simulate 
a model defined by the potential energy E[b), one introduces a set of conjugate "momenta" 
Pin and makes use of the auxiliary Hamiltonian 

HMc{p,b) = \J:pI + ^ . (10) 
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A finite-step approximation of the equations of motion arising from Hmc ^re used to guide 
the evolution of the system. A convenient choice of discretization is the leapfrog scheme 
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,N-l /x = l,2,3. (11) 



where e denotes the step size. 

The first step in the Monte Carlo procedure is to assign independent random values to the 
momenta, drawn from the distribution P[pi^) oc exp( — |p|^). Eq. 11 is then iterated n 
times, starting from the old configuration bi^. This is called one trajectory, and we denote 
the final position of the system by Pin,bj^. In the last step, the new configuration bj^ 
is accepted or rejected with probability min[l, e"^-*^^*^] for acceptance, where AHmc = 
Hmc{.P^ 1 ^^) ~ Hmc{.Pi ^)- By using the fact that the phase space map in eq. 11 is reversible 
and area-preserving, it is easily verified that this final Metropolis step makes detailed balance 
fulfilled. 

The Metropolis question involves the discretization error in the extensive quantity Hmc- 
To maintain a reasonable acceptance rate, it is therefore necessary to decrease e as iV is 
increased. The required variation of e can be estimated using fairly general arguments 
[10, 11, 12]. As the number of degrees of freedom tends to infinity, one finds that the mean 
acceptance becomes 

P„,, = (min(l,e-^^--)) = erfc ^^(A^mc)'/' j , (12) 

where {AHmc) ~ -^^^ for fixed trajectory length ne. This would imply that constant ac- 
ceptance requires e ~ N~^/^. Unfortunately, this argument fails for the models considered 
here. In fact, {AHmc): well as higher moments of AHmc-, diverges due to the singular 
nature of the energy function. As the tests below will show, the conclusion remains approx- 
imately valid in some of the cases studied, while in other a more rapid decrease of the step 
size with N is required. 

The acceptance rate depends also on the temperature, which so far was assumed to be 
fixed. This is important to note if the algorithm is used for simulated annealing. If the 
temperature is lowered, the step size must be decreased to keep the acceptance constant. 

The final choice of the two simulation parameters e and n should be based on autocorrelation 
properties. The autocorrelation function for a quantity O is given by 

Co{m) = {Omo+mOmo) " (O)' (13) 

where Omo Omo+m are measurements separated by m trajectories. The exponential 
decay expected at large m, Coim) ~ Q-'^/'^exp,o ^ defines the exponential autocorrelation 
time Teg.p O- Texp,o is the autocorrelation time of the slowest mode that couples to O. In 
our calculations, we consider instead the integrated autocorrelation time 
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which is usually much easier to measure numerically. Tint^o directly controls the statistical 
error on O , since the variance of the sample mean is given by 

m=l 

for large sample size M. 

These autocorrelation times are expressed in units of trajectories. To get a measure of 
computational effort we multiply the autocorrelation time considered, r , by n, since the CPU 
time required for each trajectory is proportional to n. Ideally, the simulation parameters 
should be chosen so as to minimize the effort E = nr for each N . This minimization requires 
extensive testing, which we have not carried out for large N . Consequently, we expect that 
our choice of parameters leads to reasonable but for large N not optimal performance. The 
corresponding effort is well described by a power law E oc , as we will see below. This 
means that the required computer time grows as JV^"*"^ , since the cost of each step in a 
trajectory scales as N"^ . 

3.2 The Gaussian Chain and Fourier Acceleration 

The dynamics of HMC can be analysed in considerable detail for the Gaussian model 
[13, 14]. We mention here a few results which motivate the Fourier acceleration method. 
For definiteness, we consider a chain with nearest-neighbour interaction only. We take the 
potential energy to be 

J =2^ ""^■^+1 = 2 ^ ^^'^ ■ ^^^^ 

i i ^1 

As before, we assume that the center of mass is held fixed and that the endpoints are free. 
The version of HMC described above is particularly easy to analyse, since E is diagonal in 
the bond variables. Using this, it is straightforward to compute the exponential autocorre- 
lation time for in the limit e ^ 0. One finds that T^xp = (— In | cosf|)~^ independent of 
i and /x, where t = ne is the trajectory length. Let us now consider HMC in the monomer 
coordinates, applied to the same model. This case can be analysed by performing the 
orthogonal transformation 

^fcM = Y]^Zl^iMCos -^^—^ k = l,...,N-l /x = 1,2,3 (17) 

1=1 

together with the same transformation of the corresponding momenta tt^^, which diagonal- 
izes 

1 1 ^"^ 1 TTk 

Hmc = 2 Zl'^^V + 2 5Z " 2 ^^^'"^ + ^fc^L) ' = 4sin2 — . (18) 

in 1=1 k II 

In the limit e ^ 0, one finds that the exponential autocorrelation time for S^.^ is given by 
Texp,k = (— In I cosa;fcf|)~^. Thus, the long-wavelength modes evolve slowly, as expected. 
The result above shows this problem can be avoided by using the bond variables. Another 
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way to speed up the evolution of the long- wavelength modes is to increase the step size for 
these, which is called the Fourier acceleration method. It is clear that this method works for 
the model in eq. 16 since the different Fourier modes evolve independently. The equations 
of motion are taken to be 



-/ — - I - _ 4 1 dE{x) 

^k^ - ^ky. + ^kT^ky. 2 T dxku , s 

k = l,...,N -1 11= 1,2,3 (19) 



^kij. — ^kii 2 T 



dEjx) , dEjx') 
Sxkii. dxkii. 



where the varying step size = e/ct'fc. Except for this change in step size the update is the 
same as in ordinary HMC. 

The Fourier acceleration technique and the use of bond variables are, in general, two distinct 
ways to improve the efficiency of HMC, and are applicable in different situations. However, 
they are equivalent for the linear chains with free endpoints considered here. In fact, 
eq. 11 can be obtained from eq. 19 by performing the same orthogonal transformation of 
itJk^kiJ. ^nd T^k,iJ.- In the calculations we have used Fourier acceleration since that makes 
the monomer coordinates more readily available. The cost of the transformation between 
monomer coordinates and Fourier variables, eq. 17, is negligble since Fast Fourier Transform 
can be utilized. 

The Fourier accelerated algorithm is set up so as to efficiently simulate the particular Gaus- 
sian model in eq. 16. The same technique can be used to adapt HMC to general Gaussian 
models. In the general case, the cost of the required coordinate transformation, corre- 
sponding to eq. 17, scales quadratically with N . In applications to models with interactions 
between all the monomer pairs, this is not necessarily severe since the cost of the energy 
computation already scales as iV^. Therefore, a study of this more general class of algo- 
rithms could be worthwhile. Especially interesting appears the possibility to make use of 
the energy function obtained in the Gaussian variational approach. We carried out some 
tests of this method for the unscreened Coulomb model, using the variational results of 
refs. [6, 7]. Somewhat disappointingly, however, these tests did not indicate any further 
gain in efficiency. 

The results for the model in eq. 16 suggest one further modification of the algorithm. In 
fact, the autocorrelation times can get large also with the Fourier accelerated algorithm, if 
t happens to be near a multiple of tt. To ensure against such accidental mode locking we 
have, following ref. [15], randomized the trajectory length t. Whether this is necessary for 
the models considered in our applications is not clear. In fact, it is possible that slightly 
better efficiency could be obtained without this randomization [16]. 



3.3 Tuning of parameters 

Next we describe our choice of simulation parameters for the models described in section 2. 
We begin with the unscreened and screened Coulomb models. Here we have randomized the 
trajectory length t = nehy drawing e from the exponential distribution with mean e, while 
holding n fixed. The average trajectory length i = ne has been taken to be independent 
of N , which for the Gaussian model would make the autocorrelation times discussed above 
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Figure 1: The acceptance probability for a) unscreened andb) screened (c^ = IM) Coulomb 
potential. The results were obtained for iV = 16 (crosses), 32 (diamonds) and 64 (squares) 
and different values of e, keeping t = 2 fixed. 



independent oi N . The remaining parameter e has been chosen so as to keep the acceptance 
constant. 

To check the behaviour of the acceptance we performed a set of calculations using iV = 16, 
32 and 64 and various e, keeping i = ne fixed. Fig. la shows the results for the unscreened 
Coulomb model plotted against e^N . Approximately, the data fall on a single curve, which 
would imply that e oc N~^/^ gives constant acceptance. This is the same scaling behaviour 
as, for example, for the Gaussian model. The results for the screened Coulomb model are 
shown in fig. lb and are different. In this case, constant acceptance requires approximately 
e oc iV~^/^. As might have been expected, the Monte Carlo dynamics seems to be more 
affected by the singularities in the potential when this is short-range. 

For the potentials with A = 2 and 2.5, we first tried to use the same tuning prescription 
with i independent of N . However, the effort E then turned out to increase faster with 
N than for the previous two models. To improve on this it was necessary to let i increase 
with iV, and we decided to take f oc N^/"^. The randomization of f has for these two models 
instead been done, for fixed e, by taking n to be uniformly distributed in 2n/3 < n < 4n/3, 
which improved the efficiency. The acceptance has been kept roughly constant by adjusting 
e, which can be done by means of short test runs. 

For all four models, we carried out a set of test runs to determine approximately optimal 
parameters for N = 32. Starting from these, we have then varied the parameters with N 
in the way described. We stress that the resulting efficiency is not expected to be optimal 
for large N. Further adjustments of f may bring significant improvements. 
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4.1(2) 




64 


5/1 


3.0 


24 


0.75 


1.44(16) 




128 


36/6 


1.2 


24 


0.75 


5.5(2) 




128 


5/1 


4.2 


42 


0.77 


1.19(11) 




256 


36/6 


1.2 
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Table 1: Details of the Monte Carlo runs for a) the unscreened and screened Coulomb 
models and b) the models with A = 2 and 2.5. "Its" indicates, in thousands, the 
total number of iterations and the number of iterations discarded for thermalization. The 
statistical error on the acceptance rate was always less than 1%. 



4 RESULTS 

In this section we present the results for the efficiency of the algorithm, which have been 
obtained using 16 < iV < 512. We also study the scaling of the end-to-end distance, and 
compare our results with the predictions mentioned in section 2. 

Details of the Monte Carlo runs are given in table 1. The longest runs required about 100 hr 
CPU time on a DEC 3000. The virial identity, see section 2, was routinely used as a check 
of the runs. The results have, when possible, been checked against those of refs. [6, 7]. The 
errors quoted on integrated autocorrelation times have been obtained by dividing the data 
into eight subsamples. Average and error for the end-to-end distance have been estimated 
through a jackknife procedure, using between 50 and 200 blocks. We checked that the errors 
were stable under change in the number of blocks. All errors given are la errors. 

4.1 Autocorrelation Times 

To monitor the efficiency of the algorithm we have mainly used the integrated autocorre- 
lation time Tint,ee for the end-to-end distance r"i^. Tint,ee gives the cost of generating one 
independent measurement of r"i^, assuming that equilibrium has been attained. In several 
cases, we also studied how the efficiency varied with the length scale considered. This was 
done by measuring the integrated autocorrelation time Tint^k for for all k. The 

results showed only a fairly weak k dependence, and the maximum Tint^k was never much 
larger than Tint,ee- No indication was found of a long autocorrelation time that is missed out 
when using Tint,ee to estimate the computational cost. We also note that the k dependence 
of Tint ,fc varies with the simulation parameters. When using shorter trajectory lengths, we 
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Figure 2: The effort E against N for a) the unscreened and screened Coulomb models and 
b) the potential with A = 2 and 2.5. 
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Table 2: Results from fits of the data for E = nTint,ee to E = cN^' . 



observed an increase in Tint^k small k. 

In fig. 2, we show the results for E = nTint,ee- The data are fairly well described by the 
power law E oc for all four models. The straight lines in the figure are fits to this form. 
The details of the fits are given in table 2. There is a clear tendency that the algorithm 
performs better for the longer range potentials, but the exponent z' does not vary much. 
The fitted values of z' all lie between 0.64 and 0.84. 



4.2 End-to-end distance 



We begin the study of the scaling with N of the end-to-end distance by forming the effective 
exponent 

1 - {rlM')l 



-In- 



N' 



(20) 



N'=2N 



21nf {rl^)N 

This gives direct information about the exponents, independently of fitting procedures. If 
the asymptotic relation in eq. 7 is valid, then we have 



V]^ K, V -\- a 



IniV 



(21) 



for large N . 
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Figure 3: The effective exponent for a) the unscreened and screened Coulomb models 
and b) for the potential with A = 2 and 2.5. 

Fig. 3 shows our results for against 1/lniV in the four different cases. For the screened 
Coulomb model the results are, although somewhat scattered, close to v = 0.6. No evidence 
is seen for corrections to asymptotic pure power-law behaviour. The data for the r~'^'^ 
potential are similar. Deviations from pure power-law behaviour are small and the values 
of are close to the predicted value v = 1j\ = 0.8. 

The situation is different in the remaining two cases. Here clear deviations from the pre- 
dicted value of V indicate the presence of corrections to the asymptotic pure power law. 
The straight line shows the predicted large N behaviour for the r"^ potential, as given by 
eq. 7 with v = \ and a = 1/2. Clearly, the agreement with this prediction is very good, 
which is confirmed by fits of the data. Using the data for 32 < iV < 512, a fit to eq. 7 yields 
V = 1.02(3) and a = 0.58(14) with P^r degree of freedom of 1. Because of the large 
statistical error on a, we also performed a fit with v = 1 fixed, which gave a = 0.498(12) 
with x^/d.f.=0.8. Hence, the results support the predicted values of both u and a. For 
the unscreened Coulomb model, we performed the same types of fits, using again data for 
32 < iV < 512. With u = I fixed we obtained a = -0.698(3) with xVd-f=3.4, while 
the fit of both the exponents gave v = 0.997(7) and a = -0.71(4) with xVd-f-=4.9. The 
quality of these fits is not perfect. The reason could be that the asymptotic form is still 
a bad approximation at iV = 32, as in fact suggested by the figure. When restricting the 
one-exponent fit to 128 < N < 512, we obtained a = —0.689(6) and indeed a very small 
X^/d.f.. The data therefore seem completely consistent with eq. 7 also for the unscreened 
Coulomb model, and suggest that a ~ —0.7. The negative sign of a implies that the average 
distance between neighbouring monomers diverges as iV ^ oo. 



5 SUMMARY 

We have developed the HMC method for simulation of single off-lattice chains with inter- 
actions between all monomer pairs. The method is exact and makes it possible to update 
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all the degrees of freedom in a computer time of order N"^ . We have tested the performance 
of the algorithm on models with short- or long-range self-repulsion. We found that these 
models can be simulated in a computer time of order JV^"*"^ with z' between 0.64 and 0.84. 
These estimates are for measurements of local as well as global quantities. The fact that 
the efficiency is similar on different length scales distinguishes HMC from currently used 
algorithms. This property makes it easier to control thermalization, which we for large N 
expect to be faster with HMC. Possible ways to further improve the efficiency of the algo- 
rithm include a systematic fine tuning of the trajectory length and the use of a higher-order 
discretization scheme [17, 18, 19]. 

We have in this paper restricted our attention to linear chains with free endpoints. We 
have discussed two versions of HMC which are equivalent for such chains, namely the 
bond variable formulation and the Fourier accelerated algorithm. The applicability to more 
general topologies is different for these two methods. Bond variables can be directly applied 
to branched structures, while the Fourier acceleration technique is well suited for the study 
of ring polymers. 

We have also presented results for the scaling of the end-to-end distance. For two of the 
models studied, we found evidence for corrections to asymptotic pure power-law behaviour. 
Our results for the potential support the asymptotic relation (r"ii\r) ~ iV'^(lniV)~" with 
V = 1 and a = 1/2, which was predicted in ref. [8]. The results for the unscreened Coulomb 
model are well described by the same expression with v = 1 and a ~ —0.7. 
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